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Abstract. Using a variational approximation we study discrete solitons of a 
nonlinear Schrodinger lattice with a cubic-quintic nonlinearity. Using an ansatz 
with six parameters we are able to approximate bifurcations of asymmetric 
solutions connecting site-centered and bond-centered solutions and resulting 
in the exchange of their stability. We show that the numerically exact and 
variational approximations are quite close for solitons of small powers. 

1. Introduction. The variational approximation (VA) has long been used as a 
semi- analytic technique to approximate solitary wave solutions of nonlinear evo- 
lution equations with an underlying Hamiltonian structure [13]. There have been 
a number of papers exploring the VA with four parameters as a relevant approx- 
imation of localized modes in discrete nonlinear Schrodinger (DNLS) equations 
[6, 14, 19]. Kaup [10] extended the variational approximation with six parameters 
that allowed him to construct not only site-centered solutions (also called on-site 
solitons) from [14] but also the bond-centered solutions (solitons centered at a mid- 
point between two adjacent sites also known as inter-site solitons). 

Site-centered and bond-centered solitons were recently considered in the context 
of the DNLS equations with competing cubic focusing and quintic defocusing non- 
linearities both in the space of one [3] and two [4] lattice dimensions. It was found 
that the two branches exchange their stability while continued with respect to the 
underlying parameters. A salient feature of this stability exchange is that the two 
branches of site-centered and bond-centered solitions do not intersect directly but 
are connected by an intermediate branch of asymmetric solitons. It was argued in 
[4] that the discrete solitons have enhanced mobility near the regimes of stability 
inversion. These properties were originally discovered in the DNLS equations with 
a saturable nonlinearity both in the space of one [8] and two [20] dimensions as 
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well as for the DNLS equations with on-site and next-site cubic nonlinearities in 
the space of one dimension [18]. It was recently shown for the same model in the 
space of two dimensions [17] that stability inversion and asymmetric solutions may 
not lead to enhanced mobility of discrete solitons if the bifurcation points of the 
solution branches are widely separated in the parameter space. 

It is the purpose of this work to apply Kaup's variational method with six pa- 
rameters from [10] to explain bifurcations of asymmetric solutions and stability 
exchange of site-centered and bond-centered discrete solitons in the context of the 
one-dimensional cubic-quintic DNLS equation. In this sense, our work is a comple- 
ment to the previous paper [3] where discrete solitons were constructed numerically 
using a dynamical reduction. Dark solitons and staggered solutions of the cubic- 
quintic DNLS equation were recently studied numerically in Refs. [15] and [16] 
respectively. 

We consider a discrete nonlinear Schrodinger equation with a cubic-quintic non- 
linearity in the form, 

#„-HC(v^„+i-hV«-i-2V'„) + B|V'nlV«-QIV'«lVn = o, nez, (1) 

where i/jn{t) : 1R.+ — > C and {C,B,Q) are real- valued parameters. Nonlinear 
Schrodinger lattices have proved to be relevant models in a variety of contexts (see 
reviews in [7, 11]), including the description of optical pulses in one-dimensional 
waveguide arrays [5]. In this application, the quantity ['0™!^ represents the intensity 
of the electric field of waveguide n, C > represents coupling strength between 
adjacent waveguides and (B, Q) measure the nonlinearity strength. A large portion 
of the literature is dedicated to Eq. (1) with (5 = 0, which would correspond to 
a medium with a Kerr nonlinearity. Recent experimental results [2, 9, 21] have 
shown that the nonlinear response of some materials is better fit with an additional 
competing quintic nonlinearity, i.e. B,Q > 0. This lends relevance to studying the 
cubic-quintic DNLS equation in the form (1). 

The cubic-quintic DNLS equation (1) has two conserved quantities, namely the 
power, 

M = ^|V„P, (2) 
nez 

and the Hamiltonian, 

nez 

Steady-state solutions have the form ipn = w„e~'^*, n e Z, where /i G M and u„ e M 
are found from the stationary DNLS equation, 

fiUn + C{un+i + Un-i ~ 2un) + Bul^ ~~ Qu^ ^ 0, n £ Z. (4) 

We seek localized solutions of the stationary DNLS equation (4) which in turn 
correspond to discrete solitons. 

Spectral stability of the steady-state solutions is studied with the linearization 
ansatz, 

V'n(t) = (un + {vn + iwn)e^* + « + iw*Je^'*^ c"'^*, n G Z, 
which leads to the spectral problem, 

1 -flWn ~ C{Wn+l + Wn-1 ~ 2Wn) ~ Bul^Wn + Quf^Wn = XVn, 
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We are looking for nonzero solutions of the linearized system in i^(Z,C^) called 
eigenvectors. The steady-state solution is called unstable if there exists at least one 
eigenvector for which Rc(A) > 0. Otherwise, the solution is called spectrally stable. 

Steady-state solutions are considered in Section 2. Spectral stability of these 
solutions is studied in Section 3. For both problems, we compare results of the 
variational approximations and the direct numerical approximations. 

2. Variational approximations of the steady-state solutions. Solutions of 
the DNLS equation (1) correspond to critical points of the Lagrangian, 

^ = E ^ - ^ndtrn)+c + - 2i^„n+|K'„i^-|ivA.r. 

nez 

(6) 

To find approximate solutions for discrete solitons, we pose a trial function in the 
form, 

V'n = ^e*'*"e-''l"-"«l, 0„ = a + fc(n-no)-f ^(n-no)', (7) 

where each of the six parameters are dependent on t. Substituting Eq. (7) into the 
Lagrangian (6) and evaluating the sums yields the effective Lagrangian, 

2 / da fcdx\ „ 42 /^dfc I3dx\^ d/3 , 



+CA^ (e'^'Sp + e-'^'Sl - 25o) + ^A^S^ - f ^^^e, (8) 



where x = 2rio — 1 and. 



coshx?7 
sinh?y ' 

cosh rj sinh XV X , 
2 sinh^ r; 2' 



Si{v,x) = — „ , 2 +0-^0, 



I ^ 1\ e X cosh ?7 sinh x?7 X% 



Si{v,x) = So{2ri,x)-, 

Sf>{-n,x) = -5*0(377, x), 



Since the center hq of the ansatz (7) can be arbitrarily chosen on [0, 1] module 
the discrete group of translations in £ Z, we shall consider x on [—1,1] (this 
restriction was used already in the derivation of (8)). The solution with x = is 
centered between lattice sites and hence is called bond- centered. The solution with 
X = ±1 is centered on a lattice site and hence is called site- centered. Solutions for 
X £ (—1, 0) U (0, 1) are called asymmetric. 

According to the variational principle, the effective Lagrangian Ccft achieves crit- 
ical values at the Euler-Lagrange equations. 



dpj dt 



cS 



0, (9) 



where pj represents a parameter of ansatz (7). Varying a yields the conservation 
law, 

A^So = M, (10) 
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Figure 1. Solutions of Eqs. (12) and (13) for C = 0.1, B = 2, and 

= 1. 



which corresponds to the dynamical invariant (2) of the power. Varying A yields, 
da dxk dkSi l3Sidx 1 d(3 S2 C 



dt 



dt 2 dt So 2 So dt 2 dt Sn 



■ e-'''S; 



2So) 



On 



(11) 



Before writing the remaining equations, it will be more convenient to make use of 
the fact that we seek steady-state solutions. This corresponds to 



f3 



^ ^ dx ^ dr; 
dt dt 







da 

and — — 

dt 



where fi is the parameter of the steady-state solution. With this assumption, the 
equations corresponding to variation of k and /3 are identically satisfied. Varying x 
and ry and making use of Eq. (11) leads to the following two equations respectively, 

A^r]smhr]x f A'^2Q coshr] coshrjx 
sinh r] cosh r/x \ 



4 cosh r; — 1 



A^ + Ce-"^ smh2T] 



0, 



(12) 



and, 



-Ce"" 1 



cosh rjx — cosh 77 x sinh rjx 



SA^Xsinh rix 



BA^ cosh 2rix 



sinh 7] coshr^X / 2 sinh 2?] cosh ryx 

4 ^'xsinhyyx , cosh 77 cosh 3?7X ^ 



(13) 



OA" 



0. 



2sinh^27; " V sinh3?7 sinh 37y 

Note these equations with = correspond to those in Ref. [10]. When % = 
(bond-centered solutions), Eq. (12) is identically satisfied and Eq. (13) becomes, 

-Ce~''(l + e-'') sinh^ 2// BA^ QA^ coshyysinh^ 2?? 



0, 



(14) 



sinh 1] 2 sinh 3?? 

which is easily solved for , giving an existence condition in terms of the parameter 
77. There exists exactly two bond-centered solutions (denoted by Bi and B2) for 
77 G (0,77cr)j which disappear as a result of the saddle-node bifurcation at jy = rjd-. 
See Fig. 1, where solutions of Eqs. (12) and (13) are plotted for C ~ 0.1, B = 2, 
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and Q = 1 with rja- « 1.56. We note, for larger C steady-state solutions become 
smoother (they approach the continuous counterpart) and so the ansatz (7), which 
is based on an exponential cusp, becomes irrelevant. 

For X 7^ Oj the term in parenthesis of Eq. (12) can be used as condition for A^, 
which can then be substituted into (13) yielding a root finding problem in (rj, x) 
parameter space. We find exactly three pairs of solutions for C = 0.1, B = 2, 
and Q — 1, which are denoted by 5*1, S'2, and Ai in Fig. 1. Branches (6*1, 6*2) and 
(^2,^1) are connected by means of the saddle-node bifurcations. Branch Si arises 
from X = as a result of the supercritical pitchfork bifurcation, while branches S2 
and Ai arise from x = as a result of the subcritical pitchfork bifurcation. 

The solution branches 5*1 and ^2 approach x = 1 rapidly, so that they correspond 
to the site-centered solitons slightly disturbed by the variational approximation. 
The solution branch Ai is a true asymmetric solution that connects the bond- 
centered and site-centered solitons. Only one branch 5*1 of site-centered solutions 
and only one branch Bi of bond-centered solutions exist in the cubic case Q = 0, 
where these two branches extend for any 77 > 0. 

To test the validity of the variational approximations we compare them against 
direct numerical solutions of the stationary DNLS equation (4). The "numerically 
exact" solutions are obtained using a Newton method with the VA solutions as 
initial seeds. The profiles of the discrete solitons that are obtained via the VA and 
the corresponding numerical solutions are shown in Fig. 2. 



If 1 




Figure 2. Plot of numerical (solid lines) and variational (markers) 
solutions for C = 0.1, B = 2, Q = 1. Top: Short bond-centered 
(left) and site-centered (right) solutions. Bottom: A taller bond- 
centered (left) and site-centered (right) solution. The asymmetric 
solution (middle) is an intermediate between the two symmetric 
profiles. The labels Ai, Bi, B2, Si, and S2 correspond to those in 
Fig. 1. 
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Figure 3. (Color online) Plot of the power M of the VA so- 
lutions (thin lines) and numerical solutions (thick lines). Stable 
branches are represented by solid blue lines and unstable branches 
by dashed red lines. There are two principle branches, one cor- 
responding to site-centered solutions (labeled Sj, j = 1..4) and 
one for bond-centered solutions (labeled Bj, j = 1..3). The two 
principle branches are connected via asymmetric solutions (labeled 
Aj,j = 1,2) at points of stability change. Only the branches 
Si, S2, Ai, Bi, and B2, are captured by the VA. Differences between 
the VA and numerical solutions are more visible in the zooms la- 
beled 'a', 'b' and 'c'. The branch labeled B3 corresponds to the 
VA that fails to capture B3. 



It is more instructive for comparison to plot the solution branches in the (/i, M) 
plane, where M is defined in Eq. (10) and fi ~ see Fig. 3. The predicted 

stability is also depicted (see Sec. 3 for details on stability computations). 

Besides the five solutions predicted within the variational approximation, there 
exists additional branches S3, S4, B3, A2 of solutions of Eq. (4) that appear on 
Fig. 3. See Fig. 4 for profiles of these "wide" solutions. The pattern of Fig. 3 
suggests existence of an infinite number of additional branches of discrete solitons 
for large power M with a characteristic snaking behavior. The snaking behavior 
was also found in the higher dimensional cubic-quintic DNLS equation [4] and the 
Swift-Hohenberg equation [1, 12]. It seems the role a higher order dispersion plays 
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Figure 4. Examples of numerical solutions shown in Fig. 3 that 
the VA did not capture. The failure of the VA is due to the wide 
nature of the solutions. Parameter values are C — 0.1, B = 2, and 
Q = 1. 



in the Swift-Hohenberg equation is replaced by discreteness in the cubic-quintic 
DNLS equation. 

The bifurcations shown in Fig. 3 look as if the asymmetric solutions connect 
exactly at the turning points of bond-centered or site-centered solutions. This is 
not the case. Rather, there exists a pair of saddle-node and pitchfork bifurcations 
and the stability exchange occurs at the pitchfork bifurcation where the asymmetric 
solutions emanate. Interestingly, the VA is able to accurately capture this subtle 
bifurcation scenario. In Fig. 5 a plot of such a bifurcation pair is shown. This 
particular pair represents a saddle- node bifurcation of a short {Bi) and tall (B2) 
bond-centered solution and a pitchfork of two asymmetric solutions (Ai) and the 
short bond-centered solution (Bi). The agreement between the variational and 
numerical approximations is impressive. The stability is also correctly predicted. 

On the other hand, bifurcations including the site-centered solutions are not cap- 
tured as well. This should not be surprising since site-centered solutions are never 
truly represented by the VA (see Fig. 1). We point out other minor discrepancies 
between the variational approximations and the numerical results in the zooms of 
Fig. 3. Zoom (a) shows where the VA solution B3 departs from the corresponding 
numerical solution B3. This is expected as the ansatz (7) is only valid for "narrow" 
solutions, i.e. those that have at most two initially excited sites. Zoom (b) shows 
that the asymmetric solution Ai is connected to the site-centered solution iS'2 (op- 
posed to 5*3) and is underestimated. The wider short site-centered solution S3 is 
not captured at all by the VA (as expected). Finally, zoom (c) shows that the VA 
falsely predicts collision of the bond-centered Bi and site-centered 6*1 solutions for 
a non-zero value of M, similar to the cubic case discussed by Kaup in [10]. 

3. Variational approximations of spectral stability. In order to determine 
stability within the VA we return to the variational equations (9) but this time 
without the assumption that P~k=^~^ = 0. Let x = (/3, x, rj, k)'^ represent 
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Figure 5. Plot of numerical (lines) and variational (markers) ap- 
proximations of power M versus /i for (7 = 0.1,-6 = 2. (5 = 1 in a 
parameter region near a pitchfork bifurcation. Stable solutions of 
the VA are indicated by circles and unstable solutions by squares. 



the four parameters of ansatz (7) after Eqs. (10) and (11) are used. To perform a 
linear stability analysis we substitute, 



X = xq + eye 



At 



A e C, 



into the four variational equations, where the steady-state solution is defined by 
a^o = (0, Xoi ?70i 0)-'" and (xc'yo) satisfy Eqs. (12) and (13). Keeping only the terms 
linear in e leads to the generalized eigenvalue problem. 



AAy = By, 

where the entries of the 4x4 matrices A and B are given by. 



an = 0, 
ai2 



Si 



ai3 = M 
ai4 = 0, 



d_ 

drj 



So 



021 = 0, 

022 = 

M d 

024 = 0, 





(Si 


d 


'S2' 




f 


KSo " 


dx 


Sq 


) 



S2 

So 



(15) 



031 

032 
033 



034 = TT 



M rdSo S2 dS 
2So [ 
0, 
0, 

M fdSoSi 
So \ dx So 



dx So dx 



dSi 
dx 



-Si 



041 = 

0.42 = 

043 = 

M 

044 — — 



M fdSoS2 _ dS2 
2So \dT] So^ dri 
0, 
0, 

M /dSoSi _ dSi 
So \ drj So di-j 
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and, 

fcl2 = 0, 622 = 0, 

613 = 0, 623 = 0, 

2CM _ K fsmhfjx ^ „ A 

bi4^—^e '(1 + S'o), ^24 = -^— — Xe '(1 + S'o) , 

00 Oo \ Sinn r] / 



&31 

632 

&33 



d_ 

dx 
d_ 
drj 



dL dSo 



dx 
dL 



dx 
dSo 



P 



dx dx 



P 



634 = 0, 



641 

&42 

^43 
644 



d_ 

dx 
d_ 
9?7 



Here we have defined, 



L = 



P = 



CM 
So 

-CM 



q2 

•^0 



{S0 



Sp — 2So) 



0. 



BM^ 



dL^dSo 
drj drj 

dL^dSo 
drj drj 



QM^ 



P 



P 



Sg ~ 2So) — 



BM' 

S'q 



^4 + 



QM^ 



Se. 



and each entry is evaluated at tlic fixed point xq. For the form of the perturba- 
tion chosen, an eigenvalue with Rc(A) > indicates instability of the steady-state 
solution a?o. 

Since A and B are real-valued matrices, both A and A* are eigenvalues. Since 
A and B has a clear block structure, if A is an eigenvalue with the eigenvector y, 
then —A is an eigenvalue with the eigenvector Sy, where S — diag(l, — 1, — 1, 1). 
Therefore, eigenvalues of (15) occur as pairs of real or imaginary eigenvalues or as 
quartets of complex eigenvalues. A typical example of eigenvalues of (15) is shown 
on the left panel of Fig. 6 for the unstable solution Bi with a pair of real and a 
pair of imaginary eigenvalues. 

The stability of the VA solutions corresponding to Fig. 3 was computed from 
eigenvalues of the generalized problem (15). The spectrum for each of the variational 
solutions shown in Fig. 3 is plotted in the left panels of Fig. 7. Only the bond- 
centered solution Bi has a branch that is predicted to be unstable. The bottom 
left panel of Fig. 7 is a zoom of small eigenvalues in the top left panel. Where the 
asymmetric solution Ai meets the bond-centered solutions Bi and B2 corresponds 
to stability exchange, whereas the connection to the site-centered solution ^2 occurs 
at A^ < and hence no stability change takes place. Note that each branch of VA 
solutions has exactly two pairs of real or imaginary eigenvalues of the generalized 
problem (15). 

A linear stability analysis of the "numerically exact" solutions was also carried 
out from the spectral stability problem (5) using the procedure described in Rcf. [3]. 
The numerical approximation of the spectrum for the unstable bond-centered solu- 
tion Bi is shown on the right panel of Fig. 6. Note that there is always a double 
zero eigenvalue in the stability problem (5) which is related to the gauge symme- 
try. This double zero eigenvalue corresponds to the eigenvector (w„,Wn) = (0,u„) 
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Figure 6. Eigenvalues corresponding to the bond-centered solu- 
tion Bi in the variational system (left) and the full system (right). 



0.4 



and the generalized eigenvector («„,?«„) = (i9^u„,0) of system (5). This double 
zero eigenvalue is captured by the variational approximation and it results in the 
conservation law (10) and the variational equation (11). On the other hand, the 
VA gives two pairs of eigenvalues and only the real pair persists in the numerical 
approximation. The purely imaginary pair of eigenvalues does not exist in the nu- 
merical solutions, where it is replaced by the continuous spectral band on the two 
symmetric intervals of the purely imaginary axis. 

Isolated, non-zero eigenvalues for solutions Bi, B2, 82, S3, and Ai are shown in 
the right panels of Fig. 7. We note that no isolated eigenvalues of Si exist. The 
overall "look" of the top panels of Fig. 7 are quite similar, but more importantly, 
the stability is correctly predicted. Note that the right panels of Fig. 7 contains 
fewer eigenvalues than the left panels of Fig. 7 because some of the purely imag- 
inary eigenvalues of the VA solutions approximate the continuous spectrum of the 
numerical solutions. An interesting difference is seen in the bottom right panel of 
Fig. 7. In the full problem, the eigenvalue pair for the asymmetric solution Ai 
vanishes when Ai intersects with the site-centered solution S'2. 

In conclusion, we have extended the results of Ref. [10] to show that not only 
does the VA faithfully represent the fundamental localized modes, but is also able 
to correctly predict the corresponding stability for small coupling constant C and 
power M. We showed this in the context of the cubic-quintic DNLS equation, 
which exhibits a family of discrete solitons, five of which were accurately captured 
by the variational approximation. It would be interesting to derive the variational 
equations in the context of time-dependent perturbations, although, the resulting 
equations would be far more complex and may undermine the utility of the method. 
It would also be relevant to extend the results of [4] and apply this analysis to the 
higher dimensional cubic-quintic DNLS equation. 
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Figure 7. (Color online) Top: plot of eigenvalues for the VA solu- 
tions (left) and numerical solutions (right) at C = 0.1, B = 2, and 
(5 = 1. Site-centered branches are plotted as dashed-dot lines (in 
various shades of red for each branch) and bond-centered branches 
as dashed lines (in various shades of blue) . The asymmetric solu- 
tion is shown as a solid black line. The thick solid black line shows 
the boundary of the continuous spectral band in the full problem. 
Bottom: zooms of the top panels near « 0. 
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